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Abstract 

The classical numerical treatment of boundary value problems defined 
on infinite intervals is to replace the boundary conditions at infinity by suit- 
able boundary conditions at a finite point, the so-called truncated boundary. 
A truncated boundary allowing for a satisfactory accuracy of the numeri- 
cal solution has to be determined by trial and errors and this seems to be 
the weakest point of the classical approach. On the other hand, the free 
boundary approach overcomes the need for a priori definition of the trun- 
cated boundary. In fact, in a free boundary formulation the unknown free 
boundary can be identified with a truncated boundary and being unknown it 
has to be found as part of the solution. 

In this paper we consider a different way to overcome the introduction of 
a truncated boundary, namely finite differences schemes defined on quasi- 
uniform grids. A quasi-uniform grid allows us to describe the infinite do- 
main by a finite number of intervals. The last node of such grid is placed 
on infinity so right boundary conditions are taken into account exactly. We 
apply the proposed approach to the Falkner-Skan model and to a problem of 
interest in foundation engineering. The obtained numerical results are found 
in good agreement with those available in literature. 
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1 Introduction 

The classical numerical treatment of boundary value problems (Bvps) on infinite 
intervals is to replace the original problem by one defined on a finite interval, 
where a finite value, the so-called truncated boundary, is used instead of infinity 
(see, for instance, Collatz [6, pp. 150-151] or Fox [fT4l p. 92]). For the accepted 
numerical solution, the value of the truncated boundary is varied until the com- 
puted results stabilize, at least, to a prefixed number of significant digits. However, 
a truncated boundary allowing for a satisfactory accuracy of the numerical solu- 
tion has to be determined by trial and errors and this seems to be the weakest point 
of the classical approach. Hence, a priori definition of the truncated boundary was 
indicated by Lentini and Keller [17] as an important area of research. 

A theory for defining asymptotic boundary conditions to be imposed at the 
truncated boundary has been developed by de Hoog and Weiss Q, Lentini and 
Keller IfTTTl and Markowich GO! I2TII . The asymptotic boundary conditions have 
been applied successfully to the numerical approximation of the so-called "con- 
necting orbits" problems of dynamical systems, see Beyn [[3l|4l[5]|. Those prob- 
lems are of interest, not only in connection with dynamical systems, but also in 
the study of traveling wave solutions of partial differential equations of parabolic 
and hyperbolic type as shown by Beyn JH, Friedman and Doedel lfT5l . Bai et al. 
El, and Liu et al. lfT9l . 

A free boundary formulation was proposed by Fazio [9] where the unknown 
free boundary was identified with a truncated boundary. In this approach the free 
boundary is unknown and has to be found as part of the solution. This free bound- 
ary approach overcomes the need for a priori definition of the truncated boundary. 
The free boundary formulation has been applied to: the Blasius problem [9], the 
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Falkner-Skan model IfTOll . a model describing the flow of an incompressible fluid 
over a slender parabola of revolution IfTTI . a connecting orbit problem lfT2ll . and a 
problem in foundation engineering [131 . 

Here we consider finite differences schemes on quasi-uniform grids applied to 
the numerical solution of BVPs defined on infinite intervals. This is a different 
way to avoid the definition of a truncated boundary. Quasi-uniform grids have 
successfully applied to the numerical solution of partial differential equations on 
unbounded domains, see Koleva lfT6ll and the references quoted therein. In the 
following sections we consider two test problems. The first is the Falkner-Skan 
model of boundary layer theory. The last one is a problem of interest in foundation 
engineering. 

We notice that the free boundary approach or a quasi-uniform grid strategy 
are as simple as the classical truncated boundary one in contrast with the asymp- 
totic boundary approach, in this context see also the opinion expressed by J. R. 
Ockendon [|23l . 

In the last section, we point out some conclusions supported by the evidences 
of the present work and indicate a possible way to extend the proposed approach 
to Bvps defined on the whole real line. 

2 Finite differences on quasi-uniform grids 

Let us consider the smooth strict monotone quasi-uniform map x = x(£), the so- 
called grid generating function, 

x=-oln(l-§), (2.1) 

where ^ E [0,1], x 6 [0,°°], and c > is a control parameter. We notice that 
*N-i — clnN for (12.11) . The problem under consideration can be discretized by 
introducing a uniform grid % n of N + 1 nodes in [0, 1] with <^o = and = ^n + h 
with h=l/N,so that x n is a quasi-uniform grid in [0, °°] . The last interval in (12.11) . 
namely [x#-i,Xjv]* is infinite but the point x N _x/ 2 is finite, because the non integer 
nodes are defined by 

x -J E - n + a \ 
x n +a — x l c, — — l , 
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with n G {0, 1, ... ,N — 1} and < (X < 1. This map allows us to describe the 
infinite domain by a finite number of intervals. The last node of such grid is 
placed on infinity so right boundary condition is taken into account correctly. 

Figured] shows the quasi-uniform mesh x = x n , n = 0, 1 , . . . , N defined by (12.11) 
with c = 5 and N = 20. 

Xq X\9 x — > oo 

oooooooo o oooooooo o o o 

5 15 

Figure 1 : A quasi-uniform mesh x = x n . We notice that the last mesh-point is 

X N = oo. 

We can define the values of u(x) on the mid-points of the grid 

x n+l~ x n+l/2 , x n+l/2~ x n 

"n+l/2 ~ u n H U n+ \ . (2.2) 

x n+l x n x n+l x n 

As far as the first derivative is concerned we can apply the following approxima- 
tion 



du 
dx 



"/i+I u n ^2 g\ 



n+l/2 2( x n+3/4- x n+l/4) 

This formulae uses the value u^ = u^,, but not x^ = °°. In order to justify finite 
difference formula (12.31) we note that, by considering u = u{£> (x)), we can write 

U n + 1 - U n 2 ( 4z+3/4 ~ £n+ 1 /4 ) ^ ^ 



du 




du 


d^ 


dx 


n-\ 


-1/2 dc, 


77+1/2 d X 



n+l/2 ^n+l—^n2(x n+3 / 4 —X n+l / 4 ) 

The last formula on the right hand side of equation (12.41) reduces to the right hand 
side of equation (12.31) because we are using a uniform grid for t, and therefore 

2 (^n+3/4 — "377+1/4) = 4W1 — 4V 

The two finite difference approximations (12.21) and (12.31) have order of accu- 
racy 0(N~ 2 ). 
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3 Bvps on infinite intervals 



Let us consider the class of Bvps defined on an infinite interval 

-^ = f(x,u) , xE [0,oo) , 

dx (3.1) 
g(u(0),uH)=0, 

where u(x) is an J— dimensional vector with ui(x) for £ = 1, . . . , d as components, 
f : [0,«>) xR J 4 R d , and giR^xR^ R d . 



A finite difference scheme on a quasi-uniform grid for the class of Bvps (13.11) 
can be defined by using the approximations given by (12.21) and (|2.3I) . We denote by 
the J— dimensional vector U„ the numerical approximation to the solution u(x n ) 



of (13.11) at the points of the mesh, that is for n = 0, 1, . . . ,N . A second order finite 



difference scheme for (13.11) can be written as follows: 

U n +i — U n — a n+ i/ 2 f (x n+ 1/2 ■,b n+ i/2^+i +c, 1+l / 2 U„) = , 

for n = 0,l,...,N-l (3.2) 
g(U ,Ujv)=0, 



where 



l n+l/2 — 2(x n+3 / 4 —x n+ i/ 4 ) 



Xn+l/2-Xn 

b n +\/2 = j ( 3 - 3 ) 

■*n+ 1 Xn 



"n+1/2 



*n+l ~~ *n+l/2 



forn = 0,l,...,iV-l. 

It is evident that (13.21) is a nonlinear system of rf- (N+ 1) equations in the d ■ 
(N+ 1) unknowns U = (Uo,Ui, . . . ,Uyv) r . We notice that b n+l / 2 ~ c n+i/2 ~ 1/2 
for all n = 0, 1, . . . ,N — 2, but when n = N—l, then b N _i/ 2 = and c^_i/ 2 = 1- 
On the contrary, we choose to set b N _x/ 2 = b N _ 3 / 2 and c N _i/ 2 = c N _i,/ 2 in order 



to avoid a suddenly jump for the coefficients of (13.21) . As it will be clear from 
the results reported in the next section this produces a much smaller error in the 
numerical solution of the system at x^. 
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For the solution of (13.21) we can apply the classical Newton's method along 
with the simple termination criterion 

where AU n g, for n = 0, 1, . . . ,N and £ = 1 , 2, . . . , d, is the difference between two 
successive iterate components and TOL is a fixed tolerance. The results listed in 
the next sections were computed by setting TOL = IE — 6. 



4 The Falkner-Skan model 

The Falkner-Skan model [8], of boundary layer theory E51 . is given by 

21 



d 3 u d 2 u 
dx 3 dx 2 

»«>) = > = o, 



i ( ^ u 

\dx 

du 







(4.1) 



dx 



(oo) = 1 , 



is a Bvp defined on a semi-infinite interval. 

As a first step we rewrite the Falkner-Skan equation in (14.11) as a first order 
system by setting 

d i u 

u i+ i{x) = —^x) , for/ = 0,1,2. 



In this way the original Bvp (14.11) specializes to 



dx 

dU2 

dx 



u 2 

«3 



— = -Ki«3-P(l -U2 ) 

Ml (0) = M2 (0)=0, M 2 H = 1 



that is, 



(4.2) 



U = (Mi,«2,W3) 
f(x,u) = («2j "3, M4, — «l«3 — P(l — U2 2 ))' 

G(u(0),uH) = ( M1 (0), M2 (0), M2 H - I) 
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4.1 Numerical results 

Some of the numerical results obtained for the test problem introduced in this 
section are reported in this subsection. Figure [2] shows the numerical solution of 
Falkner-Skan model dOT ) with P = 1. 




Figure 2: Numerical solution of Falkner-Skan model with P = 1 obtained with the 
map x = x(%) defined by d2.il) with c = 5 for N = 80. 

d^u d^u 

Table[T]lists the numerical approximations of — =-(0) and — for increas- 

dx z dx L 
ing values of N. For all values of Af we used the initial iterate 

u\(x) = U2(x) = 1/2 , ui{x) = 10 -2 . 
d 2 u 

Let us remark here that the same values of —j-^ (0) were obtained by using b N _i/ 2 = 

d 2 u 

and c N _i/ 2 = 1 in (13.21) . but, on the contrary, larger values of — ^ (°°) were com- 



puted. In table [2] we compare the obtained numerical results with those available 
in literature: the agreement is really good. 



Table 1: Numerical approximation of (0) and (oo). 



N 


iter 


d 2 u 

**<°> 


d 2 u 
dx 2 


(oo) 




20 


6 


1.238724 


-0.21 


•10" 


-7 


40 


5 


1.234124 


0.24 


•10 


-7 


80 


5 


1.232972 


-0.33 


• 10 


-7 


160 


5 


1.232684 


0.14 


•10 


-7 


320 


5 


1.232612 


-0.25 


• 10 


-7 


640 


5 


1.232594 


0.39 


• 10 


-7 


1280 


5 


1.232589 


0.33 


•10" 


-7 



Table 2: Comparison of —^(fd) and free or truncated boundary (x e and x^ re- 
spectively) for the Homann (P = 1/2) and Hiemenz (P = 1) flows. 





This paper 


Fazio una 


Nasr et al. CZ| 


Asaithambi HI 




Quasi-uniform 


Free BF 


Chebyshev method 


Finite difference 


p 


d 2 u 
x N — 2 {0) 


d 2 u 




d 2 u 




d 2 u , _ 


0.5 






3.7 


0.927805 






0.5 


oo 0.927681 


5.09 0.927680 


7.4 


0.927680 


5.67 


0.927682 


1 






3.5 


1.232617 






1 


oo 1.232589 


5.19 1.232588 


7. 


1.232588 


5.14 


1.232589 
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5 A pile in soil 



Here we consider a problem that was already used by Lentini and Keller IfTTll 
to test the asymptotic boundary conditions approach. This problem is of special 
interest here because none of the solution components is a monotone function, on 
[0,o°), see lfT2l[T3l . Let u(x) be the deflection of a semi-infinite pile embedded in 
soft soil at a distance x below the surface of the soil. The governing differential 
equation for the movement of the pile, in dimensionless form, is given by: 

d 4 u 



dx 4 



a u.„. „ a u . . du , . 

z?(0 ) = o, ^(0)=^, »H = - H = o, 

where Pi and P2 are positive material constants. As far as the boundary conditions 
are concerned, at the origin a zero moment and a positive shear P3 are assumed 
and from physical considerations it follows that u(x) and all its derivatives go to 
zero at infinity, so that, the zero asymptotic boundary conditions can be imposed. 
This problem is of interest in foundation engineering: for instance, in the design 
of drilling rigs above the ocean floor. 

The governing differential equation in (15.11) can be rewritten as a first order 
system by setting: 

d l u 

Bi + i(x)=vr(jt), for i = 0,1,2,3 . 

In this way the original Bvp (15.11) specializes to 
du\ 



u 2 , 

u 4 , (5.2) 



dx 

dU2 

dx 
du?, 

dx 

dU4 

dx 

k 3 (0)=0, u 4 (0)=P 3l kiH = 0, u 2 H = 



-Pi (l-e^ 2 " 1 ) 



that is, 



U = (wi,i<2,"3, W4) r 
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f(x, u) = (u 2 , u 3l u 4l -Pi (I- e- p ^ ) ) T 
G(u(0),uH) = (« 3 (0),«4(0) -P 3 M°°),"2H) T 



in (1341). 

5.1 Numerical results 

In this subsection we report some of the numerical results obtained for the test 
problem introduced in this section. In order to be able to compare our numerical 
results we used the same parameter values employed by Lentini and Keller IfTTll 

Pl = h Pi = \ and Pl = \- 

Moreover, we choose to consider the values of the missing initial conditions u\ (0) 
and 112(0) as representative results. 

Figure [3] shows the numerical solution of the Bvp (15.11) . 

Table [3] lists the numerical approximations of w(0) and -r-(0) for increasing 

CI 

values of N. For all values of we used the initial iterate 

u\(x) = U2(x) = u 3 (x) = «4(x) = 1 . 



Table 3: Numerical approximation ofu(0) and — (0). 

dx 



N 


iter 


u(0) 


--(o) 

dx 


20 


5 


1.420337 


0.807289 


40 


5 


1.421243 


0.807934 


80 


5 


1.421469 


0.808094 


160 


5 


1.421526 


0.808135 


320 


5 


1.421540 


0.808145 


640 


5 


1.421544 


0.80815 


1280 


5 


1.421544 


0.80815 
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Figure 3: Numerical solution of pile model d5.il) obtained with the map x = x(^) 
defined by rtZTT) with c = 5forN = 80. 



Table 4: Comparison ofu(0), and free or truncated boundary (x e and x, 

respectively) for the pile problem. 



This paper 


Fazio una 


Lentini and Keller Q/71 


Quasi-uniform = °° 


Free BFjc £ = 17.75 


Asymptotic BCs Xoo = 10 


.(0) |(0) 


.(0) |(0) 


.(0) *(0) 


1.421544 -0.808145 


1.42154 -0.808144 


1.4215 -0.80814 



As far as the missing initial conditions for the pile problem are concerned, a 
comparison of the obtained values has been considered in tabled It is easily seen 
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that our results are in good agreement with those available in literature. Lentini 
and Keller IfTTTl used the mentioned asymptotic boundary conditions and employed 
PASVAR, a routine based upon the trapezoidal difference scheme with automatic 
mesh refinement and deferred corrections as described by Lentini and Pereyra 
|[T8ll . That software is rather sophisticated because it adjusts automatically the 
mesh and the order of accuracy of the method employed. Fazio [10] used a free 
boundary formulation of the pile problem as mentioned in the introduction and 
the Keller's box finite difference method. 



6 Concluding remarks 

The numerical results for the test problems reported in the two previous sections 
show that finite difference schemes on quasi-uniform grids are an effective way 
to solve Bvps defined on infinite intervals. The application of finite difference 
schemes on quasi-uniform grids overcomes the need for a priori definition of the 
truncated boundary. 

A further advantage in using a family of quasi-uniform grids in calculations 
is the possibility to improve the obtained accuracy. The algorithm is the same for 
many finite difference methods: for numerical differentiation or integration, solv- 
ing systems of ordinary or partial differential equations. It is based on Richardson 
extrapolation, introduced by Richardson in Il24l . We carry out two calculations 
on embedded uniform or quasi-uniform grids with total number of nodes N and 
IN. All nodes of largest steps are identical to even nodes of denser grid due to 
uniformity. Suppose that we use numerical method with order of accuracy p in 
order to find an approximation of a scalar value T . For smooth enough solutions 
the error can be decomposed into a sum of inverse powers of N. The Richardson 
formula 

defines the main term of such sum. This formula is asymptotically exact when 
N — > °° if we use uniform or quasi-uniform grids. Hence, it gives the real value of 
numerical solution error without knowledge of exact solution. We can apply (16.11) 
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as a single- step correction formula 



T = T 2n + R 2N + 0(N-p- s ) 



(6.2) 



and increase the order of accuracy of our approximation. In general, we have 
s = 2 in (16.21) for a symmetrical finite difference scheme, but only s = 1 for non- 
symmetrical ones. Such enlargement of accuracy requires to solve the same prob- 
lem twice and only few further arithmetical operations and so it is very cheap. 

I could be possible to apply Richardson extrapolation to the results reported in 
table [Hand in table [3] In this way we can avoid to solve the given problem with a 
large number of grid-points. 

Let us discuss, at the end of this work, a possible way to extend the finite 
difference schemes on quasi-uniform grids to the numerical solutions of prob- 
lems defined on the whole real line, for instance, the connecting orbits problems 
mentioned in the introduction. For these problems, all boundary conditions are 
imposed at plus or minus infinity. In such a case it is possible to use the tangential 
quasi-uniform grid 



where c > is a control parameter. In fact, if n = —TV, — N+l, — 1,0, l,...,N— 
1,N, then this tangential grid cover the whole infinite line, and in particular we 
have that X-n = — °°. 
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